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Abstract. A computation scheme for solving elliptic boundary value 
problems with axially symmetric confining potentials using different sets 
of one-parameter basis functions is presented. The efficiency of the pro- 
posed symbolic-numerical algorithms implemented in Maple is shown 
by examples of spheroidal quantum dot models, for which energy spec- 
tra and eigenfunctions versus the spheroid aspect ratio were calculated 
within the conventional effective mass approximation. Critical values of 
the aspect ratio, at which the discrete spectrum of models with finite- 
wall potentials is transformed into a continuous one in strong dimen- 
sional quantization regime, were revealed using the exact and adiabatic 
classifications. 
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1 Introduction 

To analyze the geometrical, spectral and optical characteristics of quantum dots 
in the effective mass approximation and in the regime of strong dimensional 
quantization following [Ij, many methods and models were used, including the 
exactly solvable model of a spherical impermeable well [2 , the adiabatic ap- 
proximation for a lens-shaped well confined to a narrow wetting layer |3| and 
a hemispherical impermeable well [1], the model of strongly oblate or prolate 
ellipsoidal impermeable well ^ , as well as numerical solutions of the boundary 
value problems (BVPs) with separable variables in the spheroidal coordinates for 
wells with infinite and finite wall heights |6l7l8j . However, thorough comparative 
analysis of spectral characteristics of models with different potentials, including 
those with non-separable variables, remains to be a challenging problem. This 
situation stimulates the study of a wider class of model well potentials with appli- 
cation of symbolic-numerical algoritlims (SNA) and problem-oriented software, 
developed by the authors of the present paper during years |9I10I11I12I13I14| . 

Here we analyse the spectral characteristics of the following models: a spheri- 
cal quantum dot (SQD), an oblate spheroidal quantum dot (OSQD) and a prolate 
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spheroidal quantum dot (PSQD). We make use of the Kantorovich method that 
reduces the problem to a set of ordinary differential equations (ODE) [TJ]. In 
contrast to the well-known method of adiabatic representation |16] , this method 
implies neither adiabatic separation of fast and slow variables, nor the presence 
of a small parameter. We present a calculation scheme for solving elliptical BVPs 
with axially-symmetric potentials in cylindrical coordinates (CC), spherical co- 
ordinates (SC), oblate spheroidal coordinates (OSC), and prolate spheroidal co- 
ordinates (PSC). Basing on the SNA developed for axially-symmetric potentials, 
different sets of solutions are constructed for the parametric BVPs related to the 
fast subsystem, namely, the eigenvalue problem solutions (the terms and the ba- 
sis functions), depending upon the slow variable as a parameter, as well as the 
matrix elements, i.e., the integrals of the products of basis functions and their 
derivatives with respect to the parameter, which are calculated analytically by 
means of elaborated SNA MATRA, implementing in MAPLE, or numerically us- 
ing the program ODPEVP [T3] , implementing the finite-element method (FEM) . 
These terms and matrix elements form the matrices of variable coefficients in 
the set of second-order ODE with respect to the slow variable. The BVP for this 
set of ODE is solved by means of the program K ANTBP [TT] , also implementing 
the FEM. The efficiency of the calculation scheme and the SNA used is demon- 
strated by comparison of the spectra versus the ellipticity of the prolate or oblate 
spheroid in the models of quantum dots with different confining potentials, such 
as the isotropic and anisotropic harmonic oscillator, the spherical and spheroidal 
well with finite or infinite walls, approximated by smooth short-range potentials, 
as well as by constructing the adiabatic classification of the states. 

The paper is organized as follows. In Section [5] the calculation scheme for 
solving elliptic BVPs with axially-symmetric confining potentials is presented. In 
Section |3] SNA MATRA for solving parametric BVP and corresponding integrals 
implemented in Maple is described. Section 2] is devoted to the analysis of the 
spectra of quantum dot models with three types of axially-symmetric potentials, 
including the benchmark exactly solvable models. In Conclusion we summarize 
the results and discuss the future applications of our calculation scheme and the 
SNA project presented. 

2 The problem statement 

Within the effective mass approximation under the conditions of strong dimen- 
sional quantization the Schrodinger equation for the slow envelope of the wave 
function if'(f) of a charge carrier (electron e or hole h) in the models of a spher- 
ical, prolate or oblate spheroidal quantum dot (SQD, PSQD or OSQD) has the 
form 

{H - E}§{r) = {(2/ip)-ii'' + U{r) - E}W{r) = 0, (1) 

where f e is the position vector of the particle having the effective mass 
= Me (or /Xp = Hh), P — —ihVf. is the momentum operator, E is the energy 
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of the particle U{r) is the axially-symmetric potential, confining the particle 
motion in SQD. PSQD or OSQD. In Model A U{f) is chosen to be the potential 
of isotropic or anisotropic axially-symmetric harmonic oscillator with the angular 
frequency ui = 7fo^/(/^p^o)j Ifo ~ '"'^/S being an adjustable parameter: 

U^{f) = fipCj^iCiii^ + f) + C35')/2, (2) 



To = V Ci {S^o + yo) + Cs^o radius of a spherical QD (Ci = 1, C3 = 1) or that 

of a spheroidal QD (Ci = (fo/d)^, C3 = (fg/c)'*), inscribed into a spherical one, 
where d and c are the semiaxes of the ellipse which transforms into a sphere at 
d — c ~ f-Q. For Model B U (f) is the potential of a spherical or axially-symmetric 
well 

U^{r) = {0, < {i^ + f )/d^ + 52/52 < 1. [7^^ (j2 _^ y2y~2 ^ ~2/g2 > (3) 

with walls of finite or infinite height 1 <C C/q < oo- For Model C U{r) is taken 
to be a spherical or axially-symmetric diffuse potential 

U^ir) = Uo[l + exp(((i2 + f)/~a^ + j-^ - , (4) 

where s is the edge diffusiveness parameter of the function, smoothly approximat- 
ing the vertical walls of finite height C/q- Below we restrict ourselves by consider- 
ing Model B with infinite walls J/q — > oo and Model C with walls of finite height 
/Tq. We make use of the reduced atomic units: = nf? j ^^e^ is the reduced 
Bohr radius, k is the DC permittivity, E^. = Ry* = /(2^pag'^) is the reduced 
Rydberg unit of energy, and the following dimensionless quantities are intro- 
duced: #(f) = a|j-^/2^(r), 2H = H/Ry*, 2E = E/Ry*, 2U{r) = U{r)/Ry*, 
r = r/a*B, a = d/a*g, c = c/a*g, vq = fo/a*g, uj = -fro/rl = fel'/(2i?y*). For 
an electron with the reduced mass fip = fx^ — 0.067mo at k = 13.18 in GaAs: 
a% = 102A= 10.2 nm, Ry* = Er ^ 5.2 meV. 

Since the Hamiltonian H in commutes with the z-parity operator 

(z — > — z or rj — > —rj), the solutions are divided into even (ct — +1) and odd 
(a — ~1) ones. The solution of Eq. ([I}, periodical with respect to the azimuthal 
angle ip, is sought in the form of a product 'P{xf,Xs,(p) = tf'™°'(a;/, Xs)e™"''/\/27r, 
where m — 0, ±1,±2,... is the magnetic quantum number. Then the function 
(^Xf, Xg) satisfies the following equation in the two-dimensional domain f2 = 
f2,^ix,)U[2,^^ C R2\{0}, a,j{xs) = ixf^ixs),xf^^ixs)), r?., = {xf'^,xf^^): 

(Hi{xf-Xs)+H2{xs) + V{xf,Xs) - 2^;) >f'™"(x/,x,) = 0. (5) 

The Hamiltonian of the slow subsystem H2{xs) is expressed as 

Id d - 

H2{Xs) = H2{Xs) = -. — T-^g2s{xs)-^ h Vs{Xs)., (6) 

gis{Xs)axs dxs 

and the Hamiltonian of the fast subsystem Hi{xf;Xs) is expressed via the re- 
duced Hamiltonian Hf{xf] Xg) and the weighting factor gssixg)'- 

Hi{xf;xs) = g^g^{xs)Hfixf;Xs), (7) 
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Table 1. The values of conditionally fast Xf and slow Xg independent variables, 
the coefficients 5'js(a;s), gjf{xf) and the potentials Vf(xf), Vs{xs)-, Vfs{xf,Xs), in 
Eqs.©-® for SQD, OSQD and PSQD in cyhndrical (CC), spherical (SC) and 
oblate & prolate spheroidal (OSC & PSC) coordinates with (d/2)^ — ±{a'^ — c^), 
+ for OSC, - for PSC. 
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Hfi^r^^s) = - ^^gy{xf)— + Vf{xf) + Vfs{xf,Xs). 

Table [T] contains the values of conditionally fast x / and slow Xs independent 
variables, the coefficients gis(xs), g2s{xs), gssixs), 51/ (x/), g2f{xf), and the 
reduced potentials Vf{xf), Vs{xs), Vfs{xf,Xs), entering Eqs. ©-([T]) for SQD, 
OSQD and PSQD in cylindrical {x — (z, p, (p)), spherical {x = {r,ri = cos9, ip)), 
and oblate/prolate spheroidal {x — {(,,ri,(p)) coordinates [T7]. In spherical co- 
ordinates the potential V{r,ri) in Table [1] using the definitions ([5]), in the 
reduced atomic units, for Model A is expressed as 

y(r,r,) = 2r2T/(r, 77) =c^V(Ci(l -77') + Car?'), 

and for Model C as 

Vir, 77) = 2r^V{r, 77) = 2r^Uo [l + exp((r2((l - rj^)/a^ + Ca^y'/c') - l)/s)] , 

both having zero first derivatives in the vicinity of the origin r = (equlib- 
rium point). For Model B the potentials Vfs are zero, since the potential ([3]| 
is reformulated below in the form of boundary conditions with respect to the 
variables Xf and Xs- The solution of the problem ([S])-© is sought in the form 
of Kantorovich expansion [15j 

^f-(x/,x,) = X;<?r(x/;x.)x^'^'^(£;,x,), (8) 

using as a set of trial functions the eigenfunctions <?™°'(a;/; Xg) of the Hamiltonian 
Hf{xf]Xs) from ([7]), i.e., the solutions of the parametric EVP 

{Hiixf-xs) - A,(x,)} 1>r{xf-x,) = 0, (9) 
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in the interval Xf € Qxf{xs) depending on the conditionally slow variable Xg G 
as on a parameter. These solutions obey the boundary conditions 

lim Ur\^M^i) "^^'"T" '""'^ +D^r^{x.)'!>r{xf-x,\^Q m 

in the boundary points {x'J^"^ {x s) , x^^^^ {x s)} = dQxj (a^s), of the interval Qxf {xs)- 

In Eq. ((ini), Nf"'\xs) = N^f"'\ D^f"'\xs) = D'f\ unless specially declared, 

are determined by the relations n'j^"^ = 1, d'^P"^ = at to = 0, cr = +1 (or at 

(7 = 0, i.e without parity separation), A^j™*^-* — 0, Z?^""^-* = 1 at m = 0, cr = — 1 
or at m ^ 0. The eigenfunctions satisfy the orthonormality condition with the 
weighting function gif{xf) in the same interval Xf G f2xj{xs): 

{$r\<i>r) = r ^^'^^ $r(.^f;^s)<i>r(.xr,x,)g,f{xf)dxf = s,,. (n) 

Here Ai(a;s) < ... < Aj„^^(xs) < ... is the desired set of real eigenvalues. Cor- 
responding set of potential curves 2Ei{xs) < ... < 2Ej^^^{xs) < ... of Eqs. 
^ is determined by 2Ej{xs) = g:^^{xs)\j{xs). Note, that for OSC and PSC 
the desired set of real eigenvalues Xj{xs) depends on a combined parameter, 
Xs ^ = ((i/2)^2£', the product of spectral 2E and geometrical {d/2)'^ param- 
eters of the problem ([5]). The solutions of the problem (l9))- (ITT|) for Models A 
and B are calculated in the analytical form, while for Model C this is done using 
the program ODPEVP [i3j. 

Substituting the expansion into Eq. ([S]) in consideration of ^ and ([TT|) . 
we get a set of ODE for the slow subsystem with respect to the unknown vector 
functions x^"''^^\xs, E) = x^'H^s) = (xi'^Xs), ...,xflJxs)V■■ 
'^ -I-^g2s{xs):^ + 2E{xs)+IVs{xs)-2IE)x^'Hxs)^ (12) 



gis{xs) dxs dxs 
^W(x.) + ; /^"^";^Q^"-V gg4:^Q(x.)A)^W(.^). 

gis(Xs) gis{Xs) dXs 9ls{Xsj dXs J 

Here 2E(xs) — dia,g{g^^{xs)Xj{xs)), 'W{xs), and Q(a;s) are matrices of the di- 
mension jmax X jmax, 

/ N / N .d<Pi(xf:Xs)d<Pi(xf:Xs) , , „, 

W,,{xs) = W,.(xs) = y^_^^ ^ g/j ^ ' dxf, (13) 

Qij{xs) = -Qji{xs) = - / gif{xf)^,{xf;xs) — „ " dxf, 

calculated analytically for Model B and by means of the program ODPEVP [T3J 
for Model C. Note, that for Model A in SC and CS and Model B in OSC and PSC 
the variables Xf and Xg are separated, so that the matrix elements Wij{xs) = 
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Qij{xs) = are put into the r.h.s. of Eq. and Vs{xs) are substituted from 
Table [1] The discrete spectrum solutions 2E : 2Ei < 2E2 < ... < 2Et < that 
obey the boundary conditions in the points = {a;™'", x"^'^^} — dQx^ bounding 
the interval : 



lim 



(^iVi'"'^)g2.(x.)^^^^^^ + Di"^^^x^"^^PHxs)^ = 0, (14) 



where Ns"^"^^ = 1, D^™'^-' = at to = 0, cr = +1 (or at cr = 0, i.e without parity 
separation), Ns™"^ — 0, D^J^"^ — \ at m — Q^a — —1 or at m ^ 0, and the 
orthonormality conditions 

{x^'^{x,)fx^JHx,)gUxs)dx, = <5,„ (15) 

are calculated by means of the program K ANTBP [TT] . To ensure the prescribed 
accuracy of calculation of the lower part of the spectrum discussed below with 
eight significant digits we used jmax = 16 basis functions in the expansion ([S]) and 
the discrete approximation of the desired solution by Lagrange finite elements 
of the fourth order with respect to the grid pitch ^^sj-^ ^ = [^mim^k ~ -^fc-i + 

^t-i ■''max]- 



3 SNA MATRA for calculus of the BVP and integrals 

To calculate effective potentials of problem P^ - ([T5|) in each value Xg — xf, of 
the FEM grid i^^sj-^ ■) = [a;^;^, x^^xli '^^ consider a discrete representation of 
solutions (p{xf]Xs) = <?"'°' (x/ ; ) of the problem Q by means of the FEM on 
the grid, = Wo^^Lni^s), x{ = x{_^ + hi,xi= xl,^^{xs)\, in a finite 

sum: 

ftp n p 

^Xf-xs) = Y.<!>l{x,)Nl{xf) = EE<Vp(fc-i)(^^X+p(fc-i)(^/)' (16) 

^^=0 k=l r=0 

where NJ^{xf) are local functions and (!>^{xs) are node values of ^{xf^^Xs). The 
local functions NP{xf) are piece- wise polynomial of the given order p equals 
one only in the node and equals zero in all other nodes xl ^ x^ of the 
grid Sl^j^^^^y i.e., NP{x0 = 6„f^, fijV = 0,1, ... , np. The coefiicients <P^{xs) are 

formally connected with solution 'I'ixjf^.; Xg) in a node xl ~ xjf^, k — 1, . . . ,n, 
r = 0,...,p: 

The theoretical estimate for the H" norm between the exact and numerical 
solution has the order of 

\Xj{xs)-mxs)\<cih^P, \^,{xf;x,)-^';{xs)\\^<C2hP+\ (17) 
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where h-^ = niaxi<j<fi is maximum step of grid and constants ci > 0, C2 > 
do not depend on step h-^ [19\. It has been shown that we have a possibihty 
to construct schemes for solving the BVPs and integrals with high order of 
accuracy comparable with the computer one in according with the following 
estimations [131 



dx^ dx^ 



< C4hP+\ (18) 



(x,) - QlJxs)\ < c^h^P, \W,,{x,) - W!;ixs)\ < ceh-'P, (19) 



where h-^ is the grid step, p is the order of finite elements, i, j are the number 
of the corresponding solutions, and constants C3, C4, C5 and cg do not depend 
on step h-^ . Proof is straightforward following the scheme of proof of estimations 
P7|) in accordance with [19l20j . The verification of the above estimations are 
examined by numerical analysis on condensed grids and by comparison with 
examples of exact solvable models A and B. 

Let us consider the reduction of EVP ([TT|) on the interval A : x{^^^^{xs) < 
Xf < x(^^^{xs) with boundary conditions (jlOp in points x{^^^{xs) and x^^^{xs) 
rewriting in the form 

A{xs)<Pj{xf;x,) = Xj{xs)B{xs)<Pj{xf;Xs), (20) 

where A.{xs) is differential operator and B(xs) is multiplication operator are 
differentiable by parameter Xg £ ^x^- Substituting expansion (jl6p to (|20[) and 
integration with respect to Xf by parts in the interval A = yj^^^Ak, we arrive 
to a system of the linear algebraic equations 

^lAxs)1>l^{x,) = \'^{xs)hl,{x,)<Pl^{xs). (21) 

in framework of the briefly described FEM. Using p-order Lagrange elements 
[19| . we present below an Algorithm 1 for construction of algebraic problem 
(PT|) by the FEM in the form of conventional pseudocode. It MAPLE realization 
allow us show explicitly recalculation of indices /i, v and test of correspondent 
modules of parametric matrix problems, derivatives of solutions by parameter 
and calculation of integrals. 

Algorithm 1 Generation of parametric algebraic problems 

Input: 

A = yJj.^^Ak = [Mmn^^s), a^4ax(2^s)], IS interval of changing of independent vari- 
able X f , that boundaries depending on parameter Xg = xf,, ; 

f f f 
hj, = X], — x]._-^ is a grid step; 

n is a number of subintervals Ak = 
p is a order of finite elements; 

A{xs),'R{xs) are differential operators in Eq. (PO)) ; 
Output: 

NP{xf) is a basis functions in ((Tl 
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aP^^{xs)^ b^^(xs) are matrix elements in system of algebraic equations (1211) : 
Local: 

xjfj. are nodes; (f)^ ri^f) Lagrange elements; iijV — 0,1, ... ,np ; 

1: for k:=l to n do 

for r:=0 to p do 

end for; 
end for; 

3: N^^ixf) {ctfl^,{xf),Xf e A^-Q,Xf ^ A,}; 
for k:=l to n do 

for r:=l to p — 1 do 

end for; 

Nlpi^f) - e Ak;4>l+i,o{xf),Xf e Z\fe+i;0,a;/ ^ Z\fcU^fe+i}; 

end for; 

Nlp{xf) := {(j>l,p{xf),Xf £ An;0,Xf ^ Af,}; 
4: for fi, iy:=0 to do 

^U^s) := J gi{xf)NP{xf)Aixs)NEixf)dxf; 

■■= l9i{xf)NP{xf)B{xs)NP{xf)dxf; 

A 

end for; 



Remarks: 

1. For equation ^ matrix elements of the operator ([7]), and V{xf;Xs) = 
Vfs{xf,Xs) + Vf{xf) between local functions iV^(x/) and N^{xf) defined in same 
interval Aj calculated by formula, using Xf — xl_-^ + 0.5/i{(l + rjf), q,r — 0,p: 

i^i^s)),,,= 7{j^92f{Xf)i<l^iJi<t>lrY+9^^^^^^ 

+ 1 

iH^s))^^^ = J^gifixf)(t)lg(j)l,.^drif, n = q + p{k-l), v = r+p{k-l). 

2. If integrals do not calculated analytically, for example, see section |H then 
they are calculated by numerical methods [19J, by means of the Gauss quadrature 
formulae of the order p + 1 . 

3. For OSQD&PSQD model C the problem ©-([n]) has been solved using 
a grid J7j^j^^^^ [x^jj^, x^ax] = —1(20)1 (the number in parentheses denotes the 
number of finite elements of order p = 4 in each interval) . 

Generally, 10-16 iterations are required for the subspace iterations to converge 
the subspace to within the prescribe tolerance. If matrix a^ = aP{xs) in Eq. (|2ip 
is not positively defined, problem ([?!]) is replaced by the following problem: 

sP #^ = A'* bP aP = aP - ab^. (22) 



Symbolic-Numeric Algorithms for Spheroidal Quantum Dot Models 



9 



The number a (the shift of the energy spectrum) is chosen in such a way that 
matrix is positive. The eigenvector of problem (j22l) is the same, and A'' = 
a'' + a, where shift a is evaluated by the Algorithm 2. 

Algorithm 2 Evaluating the lower bound for the lowest eigenvalue of the 
generalized eigenvalue problem 

In general case it is impossible to define the lower bound for the lowest eigen- 
value ofEq. (22), because the eigenvalues Ai(a;s) < ... < ^iixs) < ■■■ < ^j^^^{xs) 
is depended on the parameter Xs- But, we can use the following algorithm to 
find the lower bound for the lowest eigenvalue Xi{xs) at fixed value of Xs- 

Step 1. Calculate LDL"^ factorization of — oB^. 

Step 2. If some elements of the diagonal matrix D are less than zero 
then put a — a ~ I and go to Step 3, else go to Step 5. 

Step 3. Calculate LDL"^ factorization of A^ — aB^. 

Step 4. If some elements of the diagonal matrix D are less than zero 
then put a = a — 1 and go to Step 3, else put a — a ~ 0.5 
and go to Step 8. 

Step 5. Put a = a + 1 and calculate L D L'^ factorization of A'' — aB^. 
Step 6. If all elements of the diagonal matrix D are greater than zero 

then put a — a + 1 and repeat Step 5. 
Step 7. Put a = a - 1.5. 
Step 8. End. 

After using the above algorithm one should find the lower bound for the 
lowest eigenvalue, and always Xi{xs) — a < 1.5 

4 Spectral characteristics of spheroidal QDs 

Models B and C for Oblate Spheroidal QD At fixed coordinate Xg of the 
slow subsystem the motion of the particle in the fast degree of freedom Xf is 
localized within the potential well having the effective width 

Lix,)^2c^l-xya^, (23) 

where L — L/a*g. The parametric BVP ©-((TT]) at fixed values of the coor- 
dinate Xs, Xs G (0,a), is solved in the interval Xf € {—L(xs)/2,L{xs)/2) 
for Model C using the program ODPEVP, and for Model B the eigenvalues 
Erio (xs) /En = 2Ei(xs), Uo = i = 1,2,..., and the corresponding parametric 
eigenfunctions (a;/; a;^), obeying the boundary conditions (jlOp and the nor- 
malization condition ifTTj) . are expressed in the analytical form: 

where the even solutions a — +1 are labelled with odd Hq = rizo + 1 = 2i — 1, and 
the odd ones cr = — 1 with even Uo = rizo + 1 = 2i, i = 1, 2, 3, ... . The effective 
potentials in Eq. for the slow subsystem are expressed analytically 
via the integrals over the fast variable Xf oi the basis functions ([M)) and their 
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a) ' b) 



Fig. 1. The energies 2E = E/Eji of even a — +1 lower states for OSQD versus 
the minor c, Cca = c/a ^ (1/5, 1) being the spheroid aspect ratio: a) well with 
impermeable walls, b) diffusion potential with 2Uq = 36, s = 0.1, the major 
semiaxis a = 2.5 and m = 0. Tine lines are minimal values 2i5J"'" = 2Ei{xs = 0) 
of potential curves. 



derivatives with respect to the parameter Xs including states with both parities 
(T = ±l: 

2E,{x,) = -^pP^ = (25) 

4c^(a^ — xj) 12 (a^ — xjy 



Wij{Xs) = 



2non:,(n2 + n;,2)(l + (-l)"<'+<) x 



2 



{rig — n'J^)'^ (a^ — x^)^ ' 

._iM!^(i±(zir:ri)^^ 



For Model B at c = a = ro the OSQD turns into SQD with known analytically 

nlm 



expressed energy levels Et = E^^ and the corresponding eigenfunctions 



^0 ?^0V^Ki+3/2(Q!n^+l,;+l/2j| 

where a„^_|_i^;+i/2 are zeros of the Bessel function of semi-integer index / + 
1/2, numbered in ascending order < an < 0:12 < ■■■ < aiy < ... by the 
integer i,v = 1,2,3,.... Otherwise one can use equivalent pairs iv -(r^ {rirj} with 
Ur = 0,1,2,... numbering the zeros of Bessel function and I = 0,1,2,... being 
the orbital quantum number that determines the parity of states a = (—1)' = 
(-1)'"(7, a = (-1)'-"* = ±1. At fixed / the energy levels Enim/Eu = 2Et, 
degenerate with respect to the magnetic quantum number m. are labelled with 
the quantum number n = n^ + l = « = 1, 2, 3, ... , in contrast to the spectrum of a 
spherical oscillator, degenerate with respect to the quantum number A = 2nr + l- 
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Figs, m [21 and [3] show the lower part of non-equidistant spectrum E{(^ca) / Eji = 
2Et and the eigenfunctions from Eq. (O for even states OSQD Models B 
and C at TO = 0. There is a one-to-one correspondence rule rio — n^o + 1 = 
2n- (1 + a)/2,ji = 1,2,3,..., rip ^ {I - \m\ - (1 - a)/2)/2, between the sets 
of spherical quantum numbers {n,l,m,a) of SQD with radius vq — a = c and 
spheroidal ones (n^ = nr,nri = I — \m\,m,a) of OSQD with the major a and 
the minor c semiaxes, and the adiabatic set of cylindrical quantum numbers 
(jizo, np, m, a) at continuous variation of the parameter (.ca — c/a. The presence 
of crossing points of the energy levels of similar parity under the symmetry 
change from spherical (ca = 1 to axial, i.e., under the variation of the parameter 
< Cca < 1, in the BVP with two variables at fixed to for Model B is caused by 
the possibility of variable separation in the OSC [T7j, i.e. the r.h.s. of Eq. (|12p 
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Fig. 3. Contour lines of the first five even-parity wave functions a — +1 in the 
xz plane of Model C of OSQD with 2Uo = 36 and s — 0.1 for the major semiaxis 
a ~ 2.5 and different values of the minor semiaxis c (^ca = c/a G (1/5, 1)). 



equals zero. The transformation of eigenfunctions, occurring in the course of a 
transition through the crossing points in Fig. [U is shown in Fig. [5] for model 
B and in Fig. [3] for model C (marked by arrows) . One can see that number 
nodes |18) the eigenfunctions ordered in according to increasing eigenvalues is 
not changed under the transition from spherical to oblate spheroidal form. So, 
at small value of deformation parameter {(ca for OSQD or (ac for PSQD) there 
are nodes only along corresponding major axis. For Model C at each value of 
the parameter a their is a finite number of discrete energy levels, limited by the 
value 2Uo of the well walls height. As shown in Fig. [TJd, the number of levels 
of OSQD, equal to that of SQD at a = c = rp, is reduced with the decrease of 
the parameter c (or Cca), in contrast to Models A and B that have countable 
spectra, and avoided crossings appear just below the threshold. 
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a) ^ b) 



Fig. 4. The energies 2E — E/Er of even cr = +1 lowest states for PSQD 
depending on the minor semiaxis a {Qac — a/c G (1/5, 1) is the spheroid aspect 
ratio): a) well with impermeable walls, b) diffusion potential, 2Uq = 36, s = 
0.1, for the major semiaxis c — 2.5 and m = Q. Tine lines are minimal values 
2£'™™ EE 2E,{xs = 0) of potential curves. 



Models B and C for Prolate Spheroidal QD In contrast to OSQD, for 
PSQD at fixed coordinate Xs of the slow subsystem the motion of the particle is 
confined to a 2D potential well with the effective variable radius 



po {x,) = a^l-xl/c\ (27) 

where po [xg) = po (xg) /as- The parametric BVP (P|)- pT|) at fixed values of 
the coordinate Xg from the interval Xg G (— c, c) is solved in the interval Xf £ 
(0,/9o {xg)) for Model C using the program ODPEVP, while for Model B the 
eigenvalues En +i {xg) jE^ = lEi {xg\ Upp + 1 = i = 1,2, and the cor- 
responding parametric basis functions (p™'^^^ {xf;Xg) = <?™(x/;Xs) without 
parity separation, obeying the boundary conditions (jlOp and the normalization 
condition (jlip . are expressed in the analytical form: 



'^|m|(^2£'„^^ + lj„| iXs)Xf) 



2 

2iJ,(^,) = %±iM, ,p™ (^^)^^_^__v • ^ (28) 



where Q:„pp+ijm| — J'^l^J^^ are positive zeros of the Bessel function of the first 
kind J|m| (x/), labeled in the ascending order with the quantum number Hpp + l = 
i = 1,2,.... The effective potentials in Ea.([T^ for the slow subsystem are 
calculated numerically in quadratures via the integrals over the fast variable x / 
of the basis functions (j28p and their derivatives with respect to the parameter Xg 
using SNA MATRA from Section 2. 

Fig. [5] illustrates the lower part of the non-equidistant spectrum E{Cac) / Eji = 
2Et of even states of PSQD Models B and C. There is a one-to-one correspon- 
dence rule Upp + 1 = Up = i = n = Ur + 1, i = 1,2, ... and Uzp = I — \m\ between 
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the sets of quantum numbers (ri, /, m, a) of SQD with the radius rp = a = c 
and spheroidal ones (n^ = 71^,11^, = I — \m\,m^a) of PSQD with the major 
c and the minor a semiaxes, and the adiabatic set of quantum numbers [n = 
Upp + l,nzp-,m, a) under the continuous variation of the parameter (^ac = a/c 
The presence of crossing points of similar-parity energy levels in Fig. Sunder the 
change of symmetry from spherical Cac = 1 to axial, i.e., under the variation of 
the parameter < Cac < 1, in the BVP with two variables at fixed m for Model 
B is caused by the possibility of variable separation in the PSC [17|, i.e. r.h.s. 
of Eq. (fT2|) equals zero. For Model C at each value of the parameter c there is 
also only a finite number of discrete energy levels, limited by the value 2Uo of 
the well walls height. As shown in Fig.[TlD the number of energy levels of PSQD, 
equal to that of SQD at a = c = ro, which is determined by the product of 
mass fie of the particle, the well depth Uq, and the square of the radius fp, is 
reduced with the decrease of the parameter d (or Cac) because of the promotion 
of the potential curve (lower bound) into the continuous spectrum, in contrast 
to Models A and B, having countable spectra. Note, that the spectrum of Model 
C for PSQD or OSQD should approach that of Model B with the growth of the 
walls height Uq of the spheroidal well. However, at critical values of the ellipsoid 
aspect ratio it is shown that in the effective mass approximation both the terms 
(lower bound) and the discrete energy eigenvalues in models of the B type move 
into the continuum. Therefore, when approaching the critical aspect ratio values, 
it is necessary to use models such as the lens-shaped self-assembled QDs with 
a quantum well confined to a narrow wetting layer [3J or if a minor semiaxis 
becomes comparable with the lattice constant to consider models (see,e.g.[5T]), 
different from the effective mass approximation. 

5 Conclusion 

By examples of the analysis of energy spectra of SQD, PSQD, and OSQD 
models with thee types of axially symmetric potentials, the efficiency of the 
developed computational scheme and SNA is demonstrated. Only Model A 
(anisotropic harmonic oscillator potential) is shown to have an equidistant spec- 
trum, while Models B and C (wells with infinite and finite walls height) possess 
non-equidistant spectra. In Model C there is a finite number of energy levels. 
This number becomes smaller as the parameter a or c (Cac or (^ca) is reduced, 
because the potential curve (lower bound) moves into the continuum. Models 
A and B have countable discrete spectra. This difference in spectra allows ver- 
ification of SQD, PSQD, and OSQD models using experimental data [2J, e.g., 
photoabsorption, from which not only the energy level spacing, but also the 
mean geometric dimensions of QD may be derived |5l7l8j . It is shown that there 
are critical values of the ellipsoid aspect ratio, at which in the approximation of 
effective mass the discrete spectrum of models with finite-wall potentials turns 
into a continuous one. Hence, using experimental data, it is possible to verify 
different QD models like the lens-shaped self-assembled QDs with a quantum 
well confined to a narrow wetting layer ^ , or to determine the validity domain 
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of the effective mass approximation, if a minor semiaxis becomes comparable 
with the lattice constant and to proceed opportunely to more adequate models, 
such as pT] , 

Note a posteriori, that the diagonal approximation of the slow- variable ODE 
(fT2|) without the diagonal matrix element Wu (so called rude adiabatic approx- 
imation) provides the lower estimate of the calculated energy levels. With this 
matrix element taken into account (adiabatic approximation) the upper estimate 
of energy is provided, unless in the domain of the energy level crossing points. 
Therefore, the Born-Oppenheimer (BO) approximation is, generally, applicable 
only for estimating the ground state at an appropriate value of the small pa- 
rameter. For Model B in the first BO approximation 2Ei w Ef^^ 4- £'|^-' is given 
by the minimal value of the slow subsystem energy E^™{xs) in the equilib- 
rium points a;^ = (namely, Ef^ = 7r^n^/(2c)^ from Eq. ([^ for OSQD and 
Ef'^ = al;^^^j^i/ from Eq. for PSQD), and by the corresponding energy 

values 2£'|"^^ — ■n{ac)~^no(2np + \m\ + l) and E} — 2(ac)^-^a„pp+i_im|(n2 + l/2) of 
the 2D and ID harmonic oscillator, respectively. In [1] it is shown that the terms 
Ei{xs) allow high-precision approximation by the Hulten potential. This can be 
accomplished by means of computer algebra software, e.g., Maple, Mathematica, 
which allows (in the rude adiabatic approximation) to obtain the lower bound of 
the spectrum by solving transcendent equations, expressed analytically in terms 
of known special functions, and to use this approach for further development of 
our SNA project. 

The software package developed is applicable to the investigation of impurity 
and exciton states in semiconductor nanostructure models. Further development 
of the method and the software package is planned for solving the quasi-2D 
and quasi-lD BVPs with both discrete and continuous spectrum, which are 
necessary for calculating the optical transition rates, channeling and transport 
characteristics in the models like quantum wells and quantum wires. 
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